// $Id: cylcreepdef.cpp 12 2008-02-13 13:00:04Z msaushkin $
/**
* @file
* Реализация функции расчёта деформаций ползучести (упругость+ползучесть) по одноосной  модели

* @author Саушкин М.Н.
*/

#include<iostream>
#include "elastblock1.hpp"
#include "creepblock1.hpp"
#include "cylcreepdef.hpp"

double cylcreepdef(double sigma, double dT, double timeMax){
   double EpsZ = 0.0; // полная осевая деформация цилиндрического образца
   double EpsT = 0.0; // полная окружная деформация цилиндрического образца
   double EpsR = 0.0; // полная радиальная деформация цилиндрического образца

   double dP   = 0.0; //приращение полной осевой деформации ползучести
   double Uk1  = 0.0; // осевая деформация ползучести uk при k=1
   double Uk2  = 0.0; // осевая деформация ползучести uk при k=2
   double dUk1 = 0.0; // приращение осевой деформации ползучести uk при k=1
   double dUk2 = 0.0; // приращение осевой деформации ползучести uk при k=2

   double Vk1  = 0.0; // осевая деформация ползучести vk при k=1
   double Vk2  = 0.0; // осевая деформация ползучести vk при k=2
   double dVk1 = 0.0; // приращение осевой деформации ползучести vk при k=1
   double dVk2 = 0.0; // приращение осевой деформации ползучести vk при k=2

   double dW  = 0.0;  // приращение осевой деформации ползучести w

   double Om  = 0.0;  // параметр повреждённости
   double dOm = 0.0;  // приращение параметра повреждённости

   double Cr  = 0.0;  // критерий повреждённости
   double dCr = 0.0;  // приращение критерия повреждённости

   double time = 0.0;        // текущее время (ч)

//   elastblock1::constinit(&elast_EI698T700); // определение парметров материала при упругости
//   creepblock1::constinit(&creep_EI698T700); // определение парметров материала при ползучести

//   FILE *stream;              //  переменная для вывода информации в файл

   // Вывод информации о задаче в файл info.txt
//   stream = freopen("info.txt", "w", stdout );
//   std::cout << "timeMax     = " << timeMax << " (ч)" << "\n";
//   std::cout << "dT          = " << dT << " (ч)" << "\n";
//   std::cout << "sigma       = " << sigma << " (kg/mm2)" << "\n";
//   elastblock1::paraminfo(); // параметры материала при упругости
//   creepblock1::paraminfo(); // параметры материала при ползучести
//   fclose(stream); // закрытие файла потока вывода

   // вычисление упругих деформаций
   EpsZ =  sigma/elastblock1::getyuong();
   EpsT = -elastblock1::getnu()*EpsZ;
   EpsR = -elastblock1::getnu()*EpsZ;
   // и запись их в файл
//   stream = freopen( "epsztr.dat", "w", stdout);
   fprintf(stdout, "time      EpsZ                     EpsT                      EpsR\n");
   fprintf(stdout, "%3.3f     %1.18f     %1.18f     %1.18f\n", time, EpsZ, EpsT, EpsR);


   do{// Вычисление деформаций ползучести в образце
      time = time+dT;
      dW   = creepblock1::calcdW(sigma, dT, Om);
      dUk1 = creepblock1::calcdUk(1, sigma, dT, Uk1, Om);
      dUk2 = creepblock1::calcdUk(2, sigma, dT, Uk2, Om);
      dVk1 = creepblock1::calcdVk(1, sigma, dT, Vk1, Om);
      dVk2 = creepblock1::calcdVk(2, sigma, dT, Vk2, Om);
      dP   = dUk1 + dUk2 + dVk1 + dVk2 + dW;

      Uk1  = Uk1 + dUk1;
      Uk2  = Uk2 + dUk2;
      Vk1  = Vk1 + dVk1;
      Vk2  = Vk2 + dVk2;

      EpsZ = EpsZ + dP;
      EpsT = EpsT - creepblock1::getmu1()*(dUk1+dUk2) - creepblock1::getmu2()*(dVk1+dVk2) - 0.5*dW;
      EpsR = EpsR - creepblock1::getmu1()*(dUk1+dUk2) - creepblock1::getmu2()*(dVk1+dVk2) - 0.5*dW;

      dOm  = creepblock1::calcdOm(sigma, dP, Om);
      Om   = Om + dOm;
      dCr  = creepblock1::calcdCr(sigma, dP, Om);
      Cr   = Cr + dCr;
      fprintf(stdout, "%3.3f     %1.18f     %1.18f     %1.18f  %1.5f\n", time, EpsZ, EpsT, EpsR, Cr);
      if (time>=timeMax) return 1; // Выход по достижению timeMax
   }
   while (Cr<=1.0);
//   fclose(stream);
   return 0; // Нормальный выход из программы
};
